!  Initial condition (density, magnetic field, velocity) 
!  for a particular configuration of a braided magnetic field
!  with magnetic domes with nulls.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: linitial_condition = .true.
!
!***************************************************************
module InitialCondition
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet  
  use Mpicomm
  use Messages
!
  implicit none
!
  include '../initial_condition.h'
!
! ampl = amplitude of the magnetic field
!
  real :: ampl = 1.0
  logical :: stretch = .false.
!
  namelist /initial_condition_pars/ &
    ampl, stretch
!
  contains
!***********************************************************************
  subroutine register_initial_condition()
!
!  Configure pre-initialised (i.e. before parameter read) variables
!  which should be know to be able to evaluate
!
!  16-december-14/simon: coded
!
!  Identify CVS/SVN version information.
!
    if (lroot) call svn_id( &
        "$Id$")
!
  endsubroutine register_initial_condition
!***********************************************************************
  subroutine initial_condition_uu(f)
!
!  Initialize the velocity field.
!
!  07-may-09/wlad: coded
!
    real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
    call keep_compiler_quiet(f)
!
  endsubroutine initial_condition_uu
!***********************************************************************
  subroutine initial_condition_lnrho(f)
!
!  Initialize logarithmic density. init_lnrho 
!  will take care of converting it to linear 
!  density if you use ldensity_nolog
!
!  07-may-09/wlad: coded
!
    real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
    call keep_compiler_quiet(f)
!
  endsubroutine initial_condition_lnrho
!***********************************************************************
  subroutine initial_condition_aa(f)
!
!  Initialize the magnetic vector potential.
!
!  16-december-14/simon: coded
!
!  Braided magnetic flux tube starting from the lower xy-plane and
!  ending at the top plane together with magnetic domes with nulls
!
!  Created 2014-12-16 by Simon Candelaresi (Iomsn)
!
    use Mpicomm, only: stop_it
    use Poisson
    use Sub
!    
    real, dimension (mx,my,mz,mfarray) :: f
!        
!   The next variables are used for the uncurling.
    integer :: l, j, ju, k
!
!   clear the magnetic field to zero
    f(:,:,:,iax:iaz) = 0.

    if (stretch .eqv. .false.) then
        do n = 1, mz, 1
            do m = 1, my, 1
                do l = 1, mx, 1
                    f(l,m,n,iax) = 0
                    f(l,m,n,iay) = ampl * (&
                        -16. * ((x(l) - 0.2) ** 2 + (y(m) - 0.2) ** 2 + (z(n) -&
                        25.8) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 0.2) **&
                        2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l)&
                        - 0.2) ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 27.6) ** 2) **&
                        (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 0.2) ** 2 + (z(n) -&
                        27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l) - 0.2) ** 2 +&
                        (y(m) + 0.2) ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.) * (x(l)&
                        - 0.2) / ((y(m) + 0.2) ** 2 + (z(n) + 25.8) ** 2) * (-z(n)&
                        - 25.8) - 18.95 * ((x(l) - 0.2) ** 2 + (y(m) + 0.2) **&
                        2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m)&
                        + 0.2) ** 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6) +&
                        18.95 * ((x(l) + 8.2) ** 2 + (y(m) - 7.8) ** 2 + (z(n) -&
                        27.6) ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m) - 7.8) **&
                        2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l)&
                        + 8.2) ** 2 + (y(m) - 8.2) ** 2 + (z(n) + 25.8) ** 2) ** (-1.&
                        / 2.) * (x(l) + 8.2) / ((y(m) - 8.2) ** 2 + (z(n) + 25.8)&
                        ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) + 8.2) ** 2 +&
                        (y(m) - 8.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)&
                        + 8.2) / ((y(m) - 8.2) ** 2 + (z(n) + 27.6) ** 2) * (-z(n)&
                        - 27.6) - 16. * ((x(l) - 7.8) ** 2 + (y(m) + 8.2) **&
                        2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m)&
                        + 8.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95&
                        * ((x(l) - 7.8) ** 2 + (y(m) + 8.2) ** 2 + (z(n) - 27.6)&
                        ** 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m) + 8.2) **&
                        2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l) -&
                        7.8) ** 2 + (y(m) + 7.8) ** 2 + (z(n) + 25.8) ** 2) ** (-1.&
                        / 2.) * (x(l) - 7.8) / ((y(m) + 7.8) ** 2 + (z(n) + 25.8)&
                        ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 7.8) ** 2 + (y(m)&
                        + 7.8) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)&
                        - 7.8) / ((y(m) + 7.8) ** 2 + (z(n) + 27.6) ** 2) * (-z(n)&
                        - 27.6) - 16. * ((x(l) - 7.8) ** 2 + (y(m) - 7.8) ** 2&
                        + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m)&
                        - 7.8) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95&
                        * ((x(l) - 7.8) ** 2 + (y(m) - 7.8) ** 2 + (z(n) - 27.6)&
                        ** 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m) - 7.8) ** 2&
                        + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + x(l) + 16. * ((x(l)&
                        - 7.8) ** 2 + (y(m) - 8.2) ** 2 + (z(n) + 25.8) ** 2) **&
                        (-1. / 2.) * (x(l) - 7.8) / ((y(m) - 8.2) ** 2 + (z(n) +&
                        25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 7.8) ** 2&
                        + (y(m) - 8.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) *&
                        (x(l) - 7.8) / ((y(m) - 8.2) ** 2 + (z(n) + 27.6) ** 2) *&
                        (-z(n) - 27.6) - 16. * ((x(l) - 0.2) ** 2 + (y(m) + 8.2)&
                        ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 0.2) /&
                        ((y(m) + 8.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8)&
                        + 18.95 * ((x(l) - 0.2) ** 2 + (y(m) + 8.2) ** 2 + (z(n) -&
                        27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) + 8.2)&
                        ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l)&
                        - 0.2) ** 2 + (y(m) + 7.8) ** 2 + (z(n) + 25.8) ** 2) **&
                        (-1. / 2.) * (x(l) - 0.2) / ((y(m) + 7.8) ** 2 + (z(n) + 25.8)&
                        ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 0.2) ** 2 +&
                        (y(m) + 7.8) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)&
                        - 0.2) / ((y(m) + 7.8) ** 2 + (z(n) + 27.6) ** 2) * (-z(n)&
                        - 27.6) - 16. * ((x(l) + 8.2) ** 2 + (y(m) + 8.2) **&
                        2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m)&
                        + 8.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) +&
                        18.95 * ((x(l) + 8.2) ** 2 + (y(m) + 8.2) ** 2 + (z(n) - 27.6)&
                        ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m) + 8.2) **&
                        2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l) +&
                        8.2) ** 2 + (y(m) + 7.8) ** 2 + (z(n) + 25.8) ** 2) ** (-1.&
                        / 2.) * (x(l) + 8.2) / ((y(m) + 7.8) ** 2 + (z(n) + 25.8)&
                        ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) + 8.2) ** 2 + (y(m)&
                        + 7.8) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)&
                        + 8.2) / ((y(m) + 7.8) ** 2 + (z(n) + 27.6) ** 2) * (-z(n)&
                        - 27.6) - 16. * ((x(l) + 8.2) ** 2 + (y(m) - 7.8) ** 2&
                        + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m)&
                        - 7.8) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) - 16.&
                        * ((x(l) + 8.2) ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 25.8)&
                        ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m) - 0.2) ** 2 +&
                        (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l) + 8.2)&
                        ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 27.6) ** 2) ** (-1.&
                        / 2.) * (x(l) + 8.2) / ((y(m) - 0.2) ** 2 + (z(n) - 27.6)&
                        ** 2) * (-z(n) + 27.6) + 16. * ((x(l) + 8.2) ** 2 + (y(m)&
                        + 0.2) ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.) * (x(l) +&
                        8.2) / ((y(m) + 0.2) ** 2 + (z(n) + 25.8) ** 2) * (-z(n)&
                        - 25.8) - 18.95 * ((x(l) + 8.2) ** 2 + (y(m) + 0.2) ** 2 +&
                        (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m)&
                        + 0.2) ** 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6) - 16.&
                        * ((x(l) - 7.8) ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 25.8) **&
                        2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m) - 0.2) ** 2 + (z(n)&
                        - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l) - 7.8)&
                        ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 27.6) ** 2) ** (-1. /&
                        2.) * (x(l) - 7.8) / ((y(m) - 0.2) ** 2 + (z(n) - 27.6) **&
                        2) * (-z(n) + 27.6) + 16. * ((x(l) - 7.8) ** 2 + (y(m) +&
                        0.2) ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.) * (x(l) - 7.8)&
                        / ((y(m) + 0.2) ** 2 + (z(n) + 25.8) ** 2) * (-z(n) -&
                        25.8) - 18.95 * ((x(l) - 7.8) ** 2 + (y(m) + 0.2) ** 2 + (z(n)&
                        + 27.6) ** 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m) +&
                        0.2) ** 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6) - 16. *&
                        ((x(l) - 0.2) ** 2 + (y(m) - 7.8) ** 2 + (z(n) - 25.8) **&
                        2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 7.8) ** 2 + (z(n)&
                        - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l) - 0.2)&
                        ** 2 + (y(m) - 7.8) ** 2 + (z(n) - 27.6) ** 2) ** (-1. / 2.)&
                        * (x(l) - 0.2) / ((y(m) - 7.8) ** 2 + (z(n) - 27.6) **&
                        2) * (-z(n) + 27.6) + 16. * ((x(l) - 0.2) ** 2 + (y(m) - 8.2)&
                        ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.) * (x(l) - 0.2)&
                        / ((y(m) - 8.2) ** 2 + (z(n) + 25.8) ** 2) * (-z(n) - 25.8)&
                        - 18.95 * ((x(l) - 0.2) ** 2 + (y(m) - 8.2) ** 2 + (z(n)&
                        + 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 8.2)&
                        ** 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6))   
                    f(l,m,n,iaz) = ampl * (&
                        -exp(-((x(l) + 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n)&
                        - 20) ** 2) / 4.) * sqrt(2.) + exp(-((x(l) - 1)&
                        ** 2) / 2. - y(m) ** 2 / 2. - ((z(n) - 12) ** 2) / 4.)&
                        * sqrt(2.) - exp(-((x(l) + 1) ** 2) / 2. - y(m) ** 2&
                        / 2. - ((z(n) - 4) ** 2) / 4.) * sqrt(2.) + exp(-((x(l)&
                        - 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n) + 4)&
                        ** 2) / 4.) * sqrt(2.) - exp(-((x(l) + 1) ** 2) / 2. -&
                        y(m) ** 2 / 2. - ((z(n) + 12) ** 2) / 4.) * sqrt(2.)&
                        + exp(-((x(l) - 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n)&
                        + 20) ** 2) / 4.) * sqrt(2.) - 16. * (((x(l)) -&
                        0.2) ** 2 + (y(m) + 8.2) ** 2 + ((z(n)) - 25.8) ** 2)&
                        ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 8.2) ** 2 +&
                        ((z(n)) - 25.8) ** 2) * (y(m) + 8.2) + 18.95 * (((x(l))&
                        - 0.2) ** 2 + (y(m) + 8.2) ** 2 + ((z(n)) - 27.6)&
                        ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 8.2)&
                        ** 2 + ((z(n)) - 27.6) ** 2) * (y(m) + 8.2) + 16. * (((x(l))&
                        - 0.2) ** 2 + (y(m) + 7.8) ** 2 + ((z(n))&
                        + 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m)&
                        + 7.8) ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 7.8) - 18.95&
                        * (((x(l)) - 0.2) ** 2 + (y(m) + 7.8) ** 2 + ((z(n))&
                        + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) /&
                        ((y(m) + 7.8) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) +&
                        7.8) - 16. * (((x(l)) + 8.2) ** 2 + (y(m) + 8.2) ** 2&
                        + ((z(n)) - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) +&
                        8.2) / ((y(m) + 8.2) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m)&
                        + 8.2) + 18.95 * (((x(l)) + 8.2) ** 2 + (y(m) + 8.2)&
                        ** 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l))&
                        + 8.2) / ((y(m) + 8.2) ** 2 + ((z(n)) - 27.6) **&
                        2) * (y(m) + 8.2) + 16. * (((x(l)) + 8.2) ** 2 + (y(m)&
                        + 7.8) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) *&
                        ((x(l)) + 8.2) / ((y(m) + 7.8) ** 2 + ((z(n)) + 25.8)&
                        ** 2) * (y(m) + 7.8) - 18.95 * (((x(l)) + 8.2) **&
                        2 + (y(m) + 7.8) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1.&
                        / 2.) * ((x(l)) + 8.2) / ((y(m) + 7.8) ** 2 + ((z(n))&
                        + 27.6) ** 2) * (y(m) + 7.8) - 16. * (((x(l)) + 8.2)&
                        ** 2 + (y(m) - 7.8) ** 2 + ((z(n)) - 25.8) ** 2) **&
                        (-1. / 2.) * ((x(l)) + 8.2) / ((y(m) - 7.8) ** 2 + ((z(n))&
                        - 25.8) ** 2) * (y(m) - 7.8) + 18.95 * (((x(l))&
                        + 8.2) ** 2 + (y(m) - 7.8) ** 2 + ((z(n)) - 27.6)&
                        ** 2) ** (-1. / 2.) * ((x(l)) + 8.2) / ((y(m) - 7.8) **&
                        2 + ((z(n)) - 27.6) ** 2) * (y(m) - 7.8) + 16. * (((x(l))&
                        + 8.2) ** 2 + (y(m) - 8.2) ** 2 + ((z(n)) +&
                        25.8) ** 2) ** (-1. / 2.) * ((x(l)) + 8.2) / ((y(m) -&
                        8.2) ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) - 8.2) - 18.95&
                        * (((x(l)) + 8.2) ** 2 + (y(m) - 8.2) ** 2 + ((z(n))&
                        + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) + 8.2) / ((y(m)&
                        - 8.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) - 8.2)&
                        - 16. * (((x(l)) + 8.2) ** 2 + (y(m) - 0.2) ** 2 +&
                        ((z(n)) - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) + 8.2)&
                        / ((y(m) - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m)&
                        - 0.2) + 18.95 * (((x(l)) + 8.2) ** 2 + (y(m) - 0.2)&
                        ** 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l))&
                        + 8.2) / ((y(m) - 0.2) ** 2 + ((z(n)) - 27.6) ** 2)&
                        * (y(m) - 0.2) + 16. * (((x(l)) + 8.2) ** 2 + (y(m)&
                        + 0.2) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l))&
                        + 8.2) / ((y(m) + 0.2) ** 2 + ((z(n)) + 25.8)&
                        ** 2) * (y(m) + 0.2) - 18.95 * (((x(l)) + 8.2) ** 2&
                        + (y(m) + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1. /&
                        2.) * ((x(l)) + 8.2) / ((y(m) + 0.2) ** 2 + ((z(n))&
                        + 27.6) ** 2) * (y(m) + 0.2) - 16. * (((x(l)) - 7.8)&
                        ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1.&
                        / 2.) * ((x(l)) - 7.8) / ((y(m) - 0.2) ** 2 + ((z(n))&
                        - 25.8) ** 2) * (y(m) - 0.2) + 18.95 * (((x(l))&
                        - 7.8) ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 27.6) **&
                        2) ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m) - 0.2) **&
                        2 + ((z(n)) - 27.6) ** 2) * (y(m) - 0.2) + 16. * (((x(l))&
                        - 7.8) ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 25.8)&
                        ** 2) ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m) + 0.2)&
                        ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 0.2) - 18.95&
                        * (((x(l)) - 7.8) ** 2 + (y(m) + 0.2) ** 2 + ((z(n))&
                        + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m)&
                        + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) + 0.2)&
                        - 16. * (((x(l)) - 0.2) ** 2 + (y(m) - 7.8) ** 2 + ((z(n))&
                        - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2)&
                        / ((y(m) - 7.8) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m)&
                        - 7.8) + 18.95 * (((x(l)) - 0.2) ** 2 + (y(m) - 7.8) **&
                        2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l))&
                        - 0.2) / ((y(m) - 7.8) ** 2 + ((z(n)) - 27.6) ** 2)&
                        * (y(m) - 7.8) + 16. * (((x(l)) - 0.2) ** 2 + (y(m) -&
                        8.2) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l))&
                        - 0.2) / ((y(m) - 8.2) ** 2 + ((z(n)) + 25.8)&
                        ** 2) * (y(m) - 8.2) - 18.95 * (((x(l)) - 0.2) ** 2 +&
                        (y(m) - 8.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1. / 2.)&
                        * ((x(l)) - 0.2) / ((y(m) - 8.2) ** 2 + ((z(n))&
                        + 27.6) ** 2) * (y(m) - 8.2) - 16. * (((x(l)) - 0.2) **&
                        2 + (y(m) - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1.&
                        / 2.) * ((x(l)) - 0.2) / ((y(m) - 0.2) ** 2 + ((z(n))&
                        - 25.8) ** 2) * (y(m) - 0.2) + 18.95 * (((x(l))&
                        - 0.2) ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 27.6) ** 2)&
                        ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) - 0.2) ** 2&
                        + ((z(n)) - 27.6) ** 2) * (y(m) - 0.2) + 16. * (((x(l))&
                        - 0.2) ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 25.8)&
                        ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 0.2)&
                        ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 0.2) - 18.95 *&
                        (((x(l)) - 0.2) ** 2 + (y(m) + 0.2) ** 2 + ((z(n))&
                        + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m)&
                        + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) + 0.2) -&
                        16. * (((x(l)) - 7.8) ** 2 + (y(m) + 8.2) ** 2 + ((z(n))&
                        - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 7.8) /&
                        ((y(m) + 8.2) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m) +&
                        8.2) + 18.95 * (((x(l)) - 7.8) ** 2 + (y(m) + 8.2) **&
                        2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l))&
                        - 7.8) / ((y(m) + 8.2) ** 2 + ((z(n)) - 27.6) ** 2) *&
                        (y(m) + 8.2) + 16. * (((x(l)) - 7.8) ** 2 + (y(m) + 7.8)&
                        ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l))&
                        - 7.8) / ((y(m) + 7.8) ** 2 + ((z(n)) + 25.8) **&
                        2) * (y(m) + 7.8) - 18.95 * (((x(l)) - 7.8) ** 2 + (y(m)&
                        + 7.8) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1. / 2.)&
                        * ((x(l)) - 7.8) / ((y(m) + 7.8) ** 2 + ((z(n)) +&
                        27.6) ** 2) * (y(m) + 7.8) - 16. * (((x(l)) - 7.8) **&
                        2 + (y(m) - 7.8) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1.&
                        / 2.) * ((x(l)) - 7.8) / ((y(m) - 7.8) ** 2 + ((z(n))&
                        - 25.8) ** 2) * (y(m) - 7.8) + 18.95 * (((x(l)) -&
                        7.8) ** 2 + (y(m) - 7.8) ** 2 + ((z(n)) - 27.6) ** 2)&
                        ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m) - 7.8) ** 2 +&
                        ((z(n)) - 27.6) ** 2) * (y(m) - 7.8) + 16. * (((x(l))&
                        - 7.8) ** 2 + (y(m) - 8.2) ** 2 + ((z(n)) + 25.8)&
                        ** 2) ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m) - 8.2) **&
                        2 + ((z(n)) + 25.8) ** 2) * (y(m) - 8.2) - 18.95 * (((x(l))&
                        - 7.8) ** 2 + (y(m) - 8.2) ** 2 + ((z(n))&
                        + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m)&
                        - 8.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) - 8.2))
                enddo          
            enddo          
        enddo
    else 
        do n = 1, mz, 1
            do m = 1, my, 1
                do l = 1, mx, 1
                    f(l,m,n,iax) = 0
                    f(l,m,n,iay) = ampl * (&
                        -16. * ((x(l) - 15.8) ** 2 + (y(m) - 15.8) ** 2 + (z(n)&
                        - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 15.8) / ((y(m) - 15.8)&
                        ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l)&
                        - 15.8) ** 2 + (y(m) - 15.8) ** 2 + (z(n) - 27.6) **&
                        2) ** (-1. / 2.) * (x(l) - 15.8) / ((y(m) - 15.8) ** 2 +&
                        (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l) - 15.8)&
                        ** 2 + (y(m) - 16.2) ** 2 + (z(n) + 25.8) ** 2) ** (-1.&
                        / 2.) * (x(l) - 15.8) / ((y(m) - 16.2) ** 2 + (z(n) + 25.8)&
                        ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 15.8) ** 2 + (y(m)&
                        - 16.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)&
                        - 15.8) / ((y(m) - 16.2) ** 2 + (z(n) + 27.6) ** 2) *&
                        (-z(n) - 27.6) + 18.95 * ((x(l) + 16.2) ** 2 + (y(m) + 16.2)&
                        ** 2 + (z(n) - 27.6) ** 2) ** (-1. / 2.) * (x(l) + 16.2)&
                        / ((y(m) + 16.2) ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6)&
                        + 16. * ((x(l) + 16.2) ** 2 + (y(m) + 15.8) ** 2 + (z(n)&
                        + 25.8) ** 2) ** (-1. / 2.) * (x(l) + 16.2) / ((y(m) +&
                        15.8) ** 2 + (z(n) + 25.8) ** 2) * (-z(n) - 25.8) - 18.95&
                        * ((x(l) + 16.2) ** 2 + (y(m) + 15.8) ** 2 + (z(n) + 27.6)&
                        ** 2) ** (-1. / 2.) * (x(l) + 16.2) / ((y(m) + 15.8) **&
                        2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6) - 16. * ((x(l) +&
                        16.2) ** 2 + (y(m) - 15.8) ** 2 + (z(n) - 25.8) ** 2) ** (-1.&
                        / 2.) * (x(l) + 16.2) / ((y(m) - 15.8) ** 2 + (z(n) -&
                        25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l) + 16.2) ** 2&
                        + (y(m) - 15.8) ** 2 + (z(n) - 27.6) ** 2) ** (-1. / 2.)&
                        * (x(l) + 16.2) / ((y(m) - 15.8) ** 2 + (z(n) - 27.6) ** 2)&
                        * (-z(n) + 27.6) + 16. * ((x(l) + 16.2) ** 2 + (y(m) - 16.2)&
                        ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.) * (x(l) + 16.2)&
                        / ((y(m) - 16.2) ** 2 + (z(n) + 25.8) ** 2) * (-z(n) -&
                        25.8) - 18.95 * ((x(l) + 16.2) ** 2 + (y(m) - 16.2) ** 2&
                        + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l) + 16.2) / ((y(m)&
                        - 16.2) ** 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6) - 16.&
                        * ((x(l) - 15.8) ** 2 + (y(m) + 16.2) ** 2 + (z(n) - 25.8)&
                        ** 2) ** (-1. / 2.) * (x(l) - 15.8) / ((y(m) + 16.2) **&
                        2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l)&
                        - 15.8) ** 2 + (y(m) + 16.2) ** 2 + (z(n) - 27.6) ** 2)&
                        ** (-1. / 2.) * (x(l) - 15.8) / ((y(m) + 16.2) ** 2 + (z(n)&
                        - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l) - 15.8) **&
                        2 + (y(m) + 15.8) ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.)&
                        * (x(l) - 15.8) / ((y(m) + 15.8) ** 2 + (z(n) + 25.8) **&
                        2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 15.8) ** 2 + (y(m)&
                        + 15.8) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)&
                        - 15.8) / ((y(m) + 15.8) ** 2 + (z(n) + 27.6) ** 2) * (-z(n)&
                        - 27.6) - 16. * ((x(l) - 15.8) ** 2 + (y(m) - 0.2) ** 2&
                        + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 15.8) / ((y(m)&
                        - 0.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95&
                        * ((x(l) - 15.8) ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 27.6)&
                        ** 2) ** (-1. / 2.) * (x(l) - 15.8) / ((y(m) - 0.2) **&
                        2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l)&
                        - 15.8) ** 2 + (y(m) + 0.2) ** 2 + (z(n) + 25.8) ** 2) **&
                        (-1. / 2.) * (x(l) - 15.8) / ((y(m) + 0.2) ** 2 + (z(n) +&
                        25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 15.8) ** 2&
                        + (y(m) + 0.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) *&
                        (x(l) - 15.8) / ((y(m) + 0.2) ** 2 + (z(n) + 27.6) ** 2)&
                        * (-z(n) - 27.6) - 16. * ((x(l) - 0.2) ** 2 + (y(m) - 15.8)&
                        ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 0.2)&
                        / ((y(m) - 15.8) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8)&
                        + 18.95 * ((x(l) - 0.2) ** 2 + (y(m) - 15.8) ** 2 + (z(n)&
                        - 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 15.8)&
                        ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l)&
                        - 0.2) ** 2 + (y(m) - 16.2) ** 2 + (z(n) + 25.8) **&
                        2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 16.2) ** 2 + (z(n)&
                        + 25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 0.2)&
                        ** 2 + (y(m) - 16.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. /&
                        2.) * (x(l) - 0.2) / ((y(m) - 16.2) ** 2 + (z(n) + 27.6)&
                        ** 2) * (-z(n) - 27.6) - 16. * ((x(l) - 0.2) ** 2 + (y(m)&
                        + 16.2) ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) -&
                        0.2) / ((y(m) + 16.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n)&
                        + 25.8) + 18.95 * ((x(l) - 0.2) ** 2 + (y(m) + 16.2) ** 2&
                        + (z(n) - 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m)&
                        + 16.2) ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16.&
                        * ((x(l) - 0.2) ** 2 + (y(m) + 15.8) ** 2 + (z(n) + 25.8)&
                        ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) + 15.8) **&
                        2 + (z(n) + 25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l)&
                        - 0.2) ** 2 + (y(m) + 15.8) ** 2 + (z(n) + 27.6) ** 2) **&
                        (-1. / 2.) * (x(l) - 0.2) / ((y(m) + 15.8) ** 2 + (z(n) +&
                        27.6) ** 2) * (-z(n) - 27.6) - 16. * ((x(l) + 16.2) ** 2 +&
                        (y(m) + 16.2) ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) *&
                        (x(l) + 16.2) / ((y(m) + 16.2) ** 2 + (z(n) - 25.8) ** 2)&
                        * (-z(n) + 25.8) - 16. * ((x(l) + 16.2) ** 2 + (y(m) - 0.2)&
                        ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) + 16.2)&
                        / ((y(m) - 0.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8)&
                        + 18.95 * ((x(l) + 16.2) ** 2 + (y(m) - 0.2) ** 2 + (z(n)&
                        - 27.6) ** 2) ** (-1. / 2.) * (x(l) + 16.2) / ((y(m) -&
                        0.2) ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l)&
                        + 16.2) ** 2 + (y(m) + 0.2) ** 2 + (z(n) + 25.8) **&
                        2) ** (-1. / 2.) * (x(l) + 16.2) / ((y(m) + 0.2) ** 2 + (z(n)&
                        + 25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) + 16.2)&
                        ** 2 + (y(m) + 0.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. /&
                        2.) * (x(l) + 16.2) / ((y(m) + 0.2) ** 2 + (z(n) + 27.6)&
                        ** 2) * (-z(n) - 27.6) - 16. * ((x(l) - 0.2) ** 2 + (y(m)&
                        - 0.2) ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) -&
                        0.2) / ((y(m) - 0.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) +&
                        25.8) + 18.95 * ((x(l) - 0.2) ** 2 + (y(m) - 0.2) ** 2 +&
                        (z(n) - 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m)&
                        - 0.2) ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. *&
                        ((x(l) - 0.2) ** 2 + (y(m) + 0.2) ** 2 + (z(n) + 25.8) **&
                        2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) + 0.2) ** 2 + (z(n)&
                        + 25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 0.2)&
                        ** 2 + (y(m) + 0.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. /&
                        2.) * (x(l) - 0.2) / ((y(m) + 0.2) ** 2 + (z(n) + 27.6) **&
                        2) * (-z(n) - 27.6) + x(l))   
                    f(l,m,n,iaz) = ampl * (&
                        -exp(-((x(l) + 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n) -&
                        20) ** 2) / 4.) * sqrt(2.) + exp(-((x(l) - 1) ** 2) / 2.&
                        - y(m) ** 2 / 2. - ((z(n) - 12) ** 2) / 4.) * sqrt(2.) - exp(-((x(l)&
                        + 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n) - 4)&
                        ** 2) / 4.) * sqrt(2.) + exp(-((x(l) - 1) ** 2) / 2. - y(m)&
                        ** 2 / 2. - ((z(n) + 4) ** 2) / 4.) * sqrt(2.) - exp(-((x(l)&
                        + 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n) + 12) ** 2)&
                        / 4.) * sqrt(2.) + exp(-((x(l) - 1) ** 2) / 2. - y(m) **&
                        2 / 2. - ((z(n) + 20) ** 2) / 4.) * sqrt(2.) - 16. * (((x(l))&
                        - 15.8) ** 2 + (y(m) - 15.8) ** 2 + ((z(n)) - 25.8) **&
                        2) ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m) - 15.8) ** 2&
                        + ((z(n)) - 25.8) ** 2) * (y(m) - 15.8) + 18.95 * (((x(l))&
                        - 15.8) ** 2 + (y(m) - 15.8) ** 2 + ((z(n)) - 27.6) ** 2)&
                        ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m) - 15.8) ** 2 + ((z(n))&
                        - 27.6) ** 2) * (y(m) - 15.8) + 16. * (((x(l)) - 15.8)&
                        ** 2 + (y(m) - 16.2) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1.&
                        / 2.) * ((x(l)) - 15.8) / ((y(m) - 16.2) ** 2 + ((z(n))&
                        + 25.8) ** 2) * (y(m) - 16.2) - 18.95 * (((x(l)) - 15.8)&
                        ** 2 + (y(m) - 16.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1.&
                        / 2.) * ((x(l)) - 15.8) / ((y(m) - 16.2) ** 2 + ((z(n)) +&
                        27.6) ** 2) * (y(m) - 16.2) - 16. * (((x(l)) + 16.2) ** 2&
                        + (y(m) - 15.8) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1. / 2.)&
                        * ((x(l)) + 16.2) / ((y(m) - 15.8) ** 2 + ((z(n)) - 25.8)&
                        ** 2) * (y(m) - 15.8) + 18.95 * (((x(l)) + 16.2) ** 2 +&
                        (y(m) - 15.8) ** 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) *&
                        ((x(l)) + 16.2) / ((y(m) - 15.8) ** 2 + ((z(n)) - 27.6) **&
                        2) * (y(m) - 15.8) + 16. * (((x(l)) + 16.2) ** 2 + (y(m)&
                        - 16.2) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l))&
                        + 16.2) / ((y(m) - 16.2) ** 2 + ((z(n)) + 25.8) ** 2)&
                        * (y(m) - 16.2) - 18.95 * (((x(l)) + 16.2) ** 2 + (y(m) -&
                        16.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1. / 2.) * ((x(l))&
                        + 16.2) / ((y(m) - 16.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m)&
                        - 16.2) - 16. * (((x(l)) - 15.8) ** 2 + (y(m) + 16.2)&
                        ** 2 + ((z(n)) - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 15.8)&
                        / ((y(m) + 16.2) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m)&
                        + 16.2) + 18.95 * (((x(l)) - 15.8) ** 2 + (y(m) + 16.2) **&
                        2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 15.8)&
                        / ((y(m) + 16.2) ** 2 + ((z(n)) - 27.6) ** 2) * (y(m) + 16.2)&
                        + 16. * (((x(l)) - 15.8) ** 2 + (y(m) + 15.8) ** 2 +&
                        ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m)&
                        + 15.8) ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 15.8)&
                        - 18.95 * (((x(l)) - 15.8) ** 2 + (y(m) + 15.8) ** 2 + ((z(n))&
                        + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m)&
                        + 15.8) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) + 15.8) - 16.&
                        * (((x(l)) - 15.8) ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) -&
                        25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m) - 0.2)&
                        ** 2 + ((z(n)) - 25.8) ** 2) * (y(m) - 0.2) + 18.95 * (((x(l))&
                        - 15.8) ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 27.6)&
                        ** 2) ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m) - 0.2) **&
                        2 + ((z(n)) - 27.6) ** 2) * (y(m) - 0.2) + 16. * (((x(l))&
                        - 15.8) ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 25.8) ** 2) **&
                        (-1. / 2.) * ((x(l)) - 15.8) / ((y(m) + 0.2) ** 2 + ((z(n))&
                        + 25.8) ** 2) * (y(m) + 0.2) - 18.95 * (((x(l)) - 15.8)&
                        ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1.&
                        / 2.) * ((x(l)) - 15.8) / ((y(m) + 0.2) ** 2 + ((z(n)) +&
                        27.6) ** 2) * (y(m) + 0.2) - 16. * (((x(l)) - 0.2) ** 2 +&
                        (y(m) - 15.8) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1. / 2.) *&
                        ((x(l)) - 0.2) / ((y(m) - 15.8) ** 2 + ((z(n)) - 25.8) **&
                        2) * (y(m) - 15.8) + 18.95 * (((x(l)) - 0.2) ** 2 + (y(m)&
                        - 15.8) ** 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l))&
                        - 0.2) / ((y(m) - 15.8) ** 2 + ((z(n)) - 27.6) ** 2) *&
                        (y(m) - 15.8) + 16. * (((x(l)) - 0.2) ** 2 + (y(m) - 16.2)&
                        ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2)&
                        / ((y(m) - 16.2) ** 2 + ((z(n)) + 25.8) ** 2) * (y(m)&
                        - 16.2) - 18.95 * (((x(l)) - 0.2) ** 2 + (y(m) - 16.2) **&
                        2 + ((z(n)) + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) /&
                        ((y(m) - 16.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) - 16.2)&
                        - 16. * (((x(l)) - 0.2) ** 2 + (y(m) + 16.2) ** 2 + ((z(n))&
                        - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m)&
                        + 16.2) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m) + 16.2) + 18.95&
                        * (((x(l)) - 0.2) ** 2 + (y(m) + 16.2) ** 2 + ((z(n))&
                        - 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 16.2)&
                        ** 2 + ((z(n)) - 27.6) ** 2) * (y(m) + 16.2) + 16. * (((x(l))&
                        - 0.2) ** 2 + (y(m) + 15.8) ** 2 + ((z(n)) + 25.8)&
                        ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 15.8) **&
                        2 + ((z(n)) + 25.8) ** 2) * (y(m) + 15.8) - 18.95 * (((x(l))&
                        - 0.2) ** 2 + (y(m) + 15.8) ** 2 + ((z(n)) + 27.6) ** 2)&
                        ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 15.8) ** 2 + ((z(n))&
                        + 27.6) ** 2) * (y(m) + 15.8) - 16. * (((x(l)) + 16.2)&
                        ** 2 + (y(m) + 16.2) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1.&
                        / 2.) * ((x(l)) + 16.2) / ((y(m) + 16.2) ** 2 + ((z(n))&
                        - 25.8) ** 2) * (y(m) + 16.2) + 18.95 * (((x(l)) + 16.2)&
                        ** 2 + (y(m) + 16.2) ** 2 + ((z(n)) - 27.6) ** 2) ** (-1.&
                        / 2.) * ((x(l)) + 16.2) / ((y(m) + 16.2) ** 2 + ((z(n)) -&
                        27.6) ** 2) * (y(m) + 16.2) + 16. * (((x(l)) + 16.2) ** 2&
                        + (y(m) + 15.8) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.)&
                        * ((x(l)) + 16.2) / ((y(m) + 15.8) ** 2 + ((z(n)) + 25.8)&
                        ** 2) * (y(m) + 15.8) - 18.95 * (((x(l)) + 16.2) ** 2 +&
                        (y(m) + 15.8) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1. / 2.) *&
                        ((x(l)) + 16.2) / ((y(m) + 15.8) ** 2 + ((z(n)) + 27.6) **&
                        2) * (y(m) + 15.8) - 16. * (((x(l)) + 16.2) ** 2 + (y(m)&
                        - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1. / 2.) * ((x(l))&
                        + 16.2) / ((y(m) - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) *&
                        (y(m) - 0.2) + 18.95 * (((x(l)) + 16.2) ** 2 + (y(m) - 0.2)&
                        ** 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l)) + 16.2)&
                        / ((y(m) - 0.2) ** 2 + ((z(n)) - 27.6) ** 2) * (y(m)&
                        - 0.2) + 16. * (((x(l)) + 16.2) ** 2 + (y(m) + 0.2) ** 2 +&
                        ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l)) + 16.2) / ((y(m)&
                        + 0.2) ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 0.2) -&
                        18.95 * (((x(l)) + 16.2) ** 2 + (y(m) + 0.2) ** 2 + ((z(n))&
                        + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) + 16.2) / ((y(m)&
                        + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) + 0.2) + 18.95&
                        * (((x(l)) - 0.2) ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 27.6)&
                        ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) - 0.2) **&
                        2 + ((z(n)) - 27.6) ** 2) * (y(m) - 0.2) + 16. * (((x(l))&
                        - 0.2) ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 25.8) ** 2)&
                        ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 0.2) ** 2 + ((z(n))&
                        + 25.8) ** 2) * (y(m) + 0.2) - 18.95 * (((x(l)) - 0.2)&
                        ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1.&
                        / 2.) * ((x(l)) - 0.2) / ((y(m) + 0.2) ** 2 + ((z(n)) + 27.6)&
                        ** 2) * (y(m) + 0.2) - 16. * (((x(l)) - 0.2) ** 2 + (y(m)&
                        - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1. / 2.) * ((x(l))&
                        - 0.2) / ((y(m) - 0.2) ** 2 + ((z(n)) - 25.8) ** 2)&
                        * (y(m) - 0.2))
                enddo          
            enddo          
        enddo
    end if
!
  endsubroutine initial_condition_aa
!***********************************************************************
    subroutine read_initial_condition_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=initial_condition_pars, IOSTAT=iostat)
!
    endsubroutine read_initial_condition_pars
!***********************************************************************
    subroutine write_initial_condition_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=initial_condition_pars)
!
    endsubroutine write_initial_condition_pars
!***********************************************************************
!********************************************************************
!************        DO NOT DELETE THE FOLLOWING       **************
!********************************************************************
!**  This is an automatically generated include file that creates  **
!**  copies dummy routines from noinitial_condition.f90 for any    **
!**  InitialCondition routines not implemented in this file        **
!**                                                                **
  include '../initial_condition_dummies.inc'
!********************************************************************
endmodule InitialCondition
